Nonstationary Velocity Estimation 

By T. M. BURFORD 

(Manuscript received January 3, 1958) 

A nonstationary noise may frequently be approximated by the product 
of a stationary noise and a deterministic function of time. From observa- 
tions of tlie sum of such a nonstationary noise and a linear signal, an esti- 
mate of the rate of change of the signal is found . More exactly, a random 
function, x(t), is assumed to be one of the following : 

a + bt + g(t)n(t), 

a + bt + g(t) I"' h(t - r)n(r) dr, or 

J— 00 

a + bt + I h(t - r)g(r) ll (r) dr, 

J— 00 

where a and b are constants,h(t) is the impulse response of a lumped parame- 
ter filter, n (t) is white noise and g(t) is a nonzero deterministic function. 
A least squares estimate of b is found as a linear operation on a finite sam- 
ple of x(t). 

I. WHITE NOISE CASE 

A well-known estimation problem is that of forming a mean square 
estimate of the constant, b, in the random function x(t) = a + bt + n{t), 
where /?.(*) is white noise and a and b are unknown constants. The esti- 
mate of b is required to be the linear operation 



b = f K(t - r).r(r) dr, 

Jt — T 



where K(z) vanishes outside (0, T) and b is to equal b in the absence of 
noise. The solution is known ' to be 

K{z) = ^ (T - 2f); £ s £ T 

= 0: elsewhere. 

1009 
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The above problem will be generalized here to include one type of non- 
stationary noise. 
Let 

x(t) = a + bt + g(t)n(t), (2) 

where a and b are unknown constants,* n{t) is white noise of unit spectral 
density and g(t) is a nonzero, deterministic function. Formally, it fol- 
lows that 

E\g{t)n(t)g{t')n{i')] = g(()g(l')8(t - t'). (3) 

We wish an estimate of 6 in the form 



b = [ K(t, t - t)x(t) (It, 

Jt-T 



(4) 



where K(t, z) vanishes for z outside (0, 7'). A constraint is that, in the 
absence of noise, (4) should give b exactly. Therefore 

b = J K(t, I - t)((i + In) <h, 

which implies, with a change of variable, that 

( K(t, z) dz = 0, (5) 

Jo 

and 

zK(t,z) dz = -1. (G) 



L 



Using (3), (5) and (6) and again changing variables we find that the ex- 
pected error is now 

E(b -b) 2 = f K 2 (t, z)g\t - z) dz. (7) 

Jo 

The minimization of (7) is easily done by the usual variational tech- 
nique, using (5) and (6) as isoperimetric constraints. The variation gives 

2g\t - z)K( t, z) + X + \iz = 0. 

Therefore 

* The linear a + bt is used here for simplicity. The coefficients of a general 
polynomial may be estimated in a similar way. 
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where X unci m are undetermined multipliers to be fixed by (5) and (G). 
From (5), 

[ T „, v , X f T dz m f T zdz n 

/ K(t, z) dz = -- - 1Tt r - £ / -rr. r = 0, 

Jo ' 2 Jo (J-{t - z) 2 Jo g 2 (t - z) 

and, from (6), 

C v(* \ , X f T zdz " f r ±J* 1 

/ zK(t, z) dz = — - / - T7 - l — r — x I -TJ-, x = — 1. 

Jo 2 Jo g 2 (t, z) 2 Jo g-(t - z) 

Solving the preceding equations for X and n and substituting into (8) 
gives 

K^' M-Mi- IY °*'* T (9) 

= 0; elsewhere, 

with 

... r z j dz 

•'o (/-{t — z) 
which gives the minimized error 

E «> ~ ^ = A A A ° A* ' 

A A 2 - Af- 

The denominator of K(t, z) does not vanish, since g(t — z) is assumed 
nonzero, and 

( r dz>0 i 




A A 2 - Ai 2 > 0. 

In general, the K(t, z) found here defines a time variable filter over a 
finite time interval. However, if g is a constant, (9) reduces to (1). An- 
other interesting special case is that in which g is assumed to be of ex- 
ponential form, which implies that 

g{t - z) = g(t)g(-z), 
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and therefore 
Substitution into (9) gives 



= 0; elsewhere. 

Now, however, the B/s are constants instead of functions of t, as the 
Aj's were. For an exponential g, therefore, the filter given by (9) is time- 
invariant. 

The function g{t) is not necessarily continuous. For example, let 

g(x) = g\\ x < fi 

= 02 ; x > fi. 



Then 



1 T i+l 



1 T 



1+1 



■j + 1' 



> t, 



which implies that K(t, z) is linear but has a discontinuity in value and 
slope at z = fi. 

Many applications involve a linear g(t). For this case, it is possible to 
plot a one-parameter family of curves which completely describe K(t, z). 
If g(t) equals a + fit, we define Q = — (a + fit) /fit. Then, by a change 
of variable, 

A ' " ? So (QTi)" 2 " /3 2 Gj ' 

and substitution into (9) gives 

i**fc-)- r /at 
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All of the time dependence is now included in Q. Fig. 1 shows T'K(t, z) 
plotted against z/T for several values of Q. A positive value of Q implies 
that | g(t) | is decreasing with time and therefore that recent data are 
superior. Values of Q in (-1, 0) are not permitted, since such a Q would 
imply that g(l) vanished somewhere in (t - T,t). If g(t) is constant, Q 
is infinite and (1) applies. 
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Fig. 1 — Weighting functions, K{t, z), for linear g{t). 
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The K(t, z) just found for the linear g{t) case may be substituted into 
(7) to give the variance of b as a function of Q : 



| E(b - b) 2 = 



t-(log£±iYQ(Q + l)]\ 



(10) 



An interesting comparison may be made with the variance resulting 
from using K(z) of (1) in (7), which gives 

!#(6-&) a = 3(2Q + l) 2 + ?. (11) 

0- 5 

Equation (11) shows the error to be expected from using a K(z), 
which is optimum in the stationary case, on this type of nonstationary 
data. 

Equations (10) and (11) are plotted in Fig. 2. The difference between 
the curves is a measure of the improvement possible* with time-variable 
smoothing. 

As Q approaches <x> , (10) may be written 



T - E(b - b) 2 = 3(2Q + l) 2 - '- + 



I- + (-o> 



therefore (10) and (11) will asymptotically differ by 16/5. 

Fig. 2 is plotted for Q > 0; however, it may be used for Q < — 1 since 
both (10) and (11) are even on either side of Q = —0.50. 

Fig. 1 also suggests a simple approximate realization for K((, z). The 
several curves have approximately common intersections near z/T = 0.2 
and 0.8. Therefore, K{t, z) might be represented as 

K(t,z) = K x {z) + f(Q)K 2 (z), 

where Ki(z) is one of the K(t, z) curves near the middle of the range of 
Q's to be considered. The function K 2 (z) vanishes near the points 
z/T = 0.2 and 0.8 and takes on relatively small values elsewhere. This 
approximate realization thus involves two time-invariant filters and a 
multiplier which depends on Q and thereby on t. 

II. PRE-FILTERING OF NOISE 

Filtering and multiplication by g(t) are not, in general, commutative 
operations, so a choice must be made between 

0(0 / h(t - r)ft(r) dr (12) 

J— on 



* It has been shown by E. N. Gilbert that, even if a linear g(t) vanishes in 
the smoothing interval, the best possible estimate has variance P 2 /T and, there- 
fore, perfect estimation is not possible. 
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Fig. 2 — Variance for linear g(t). 
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and 



I h(t - T)g(r)n(T) dr 

J— QO 



(13) 



us noise functions to be added to a + bt. In certain cases, the representa- 
tions (12) and (13) are virtually equivalent in the sense that, given 
hit - t), there exists an h(l - r) such that 

g(t) f h{t - t)?i(t) dr = f h(t - r)g{r)n{r) dr. 

From the point-wise uneorrelation of n(r), the preceding equation, if it 
is true at all, must be true for all values of t. Therefore 

K(t - r) = ^ Ht - r), 

which implies that g{t)/g(r) must be a function of (/ - r). The only 
nonconstant g(t) satisfying this condition is the exponential. If git) is an 
exponential, the choice between (12) and (13) is arbitrary. In general, 
however, the choice is not arbitrary and depends on physical considera- 
tions. The forms (12) and (13) will now be treated as Cases 1 and 2. 
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Case 1 

The observed quantity, x(l), is 

x(t) = a + bt + g(t) I h(t - t)u{t) dr, 

j J— 00 

where u{t) is white noise of spectral density unity. The transform of 
h{t) is assumed to be 

m 

T[h(t)] = -± ; m < n. 



It is also assumed that the numerator and denominator of T[h(t)\ have 
neither common factors nor multiple roots. The covariance of 

J h(t - t>(t) dr 

is denned to be p(r). A particular K(t, z) is sought such that 

b = I K(t, t - t).v(t) dr 

Jt—T 

minimizes 

E(b -b) 2 = f dz g(t - z)K(t, z) I dz' g(t - z')K(t, z') P {z - /). (14) 
•»o Jo 

In obtaining (14), it was assumed that, in the absence of noise, b 
equals b, or 

f K(t,z) dz = [ zK(l,z) dz = -1. (15) 

Jo Jo 

Using the constraints of (15), the variation of (14) gives the integral 
equation* 

f dz'K{t, z')g{t - z') P {z - z') = A+_^ . (16) 

Jo g(t - z) 

From the form of T[h(t)], it is evident that 

p{z ~ z) = ^l„A^o e ■ du > 



* A similar equation, of somewhat different origin, has been considered bv 
Ule. 3 
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where A/(w 2 ) and N(o?) are polynomials in a of degrees m and n. Now 
form 4 the operation N(-d 2 /dz 2 ) and apply it to (16), giving 

h f *«* ^ - /] L *<*•"" *• ■ N (-*) w^> ■ 

which is, formally, the same as 
or, the differential equation 

*(-£)**#-*-*(-£)#$• (17) 

The general bounded solution of (17) will not necessarily satisfy (16), 
however. As is well known, 2 singularity functions up to order n - m - 1 
should be added at the end points of the interval in order to satisfy the 
boundary conditions of the integral equation (16). 
Therefore, let 

KIM - *) = ?-""- +/> " H-QU^y (ig) 

+ "if M°'( 2 ) + l" c/ j) (z - T), 

II 

where 5 <0) is the Dirac delta function, 5 (1> its derivative, etc. The func- 
tion k(z - x) is the Green's function associated with M(-(d 2 /dz 2 )) and 
the 2m constants m,- are the roots of il/(-X 2 ). The a, , bj and cj com- 
prise 2n undetermined constants. They may be determined by substitut- 
ing (18) into the integral equation (16) and requiring an identity in the 
2n exponentials of p. For example, assume h(x) = e ", which implies 
that 

M(u) = _J 

N(u-) a 2 + or ' 

Then 

Kit, z)g(t - .) - (a ! - g £±J| + WW + crf( Z - f>, 
and substitution into (16) gives the following equations for b and c : 
/ d\ X + M g c a -(a + -\ X + MZ 
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Except for evaluating X and n from (15), the solution is 
K(t, z) = 

Case 2 

The observed quantity x(t) now is 

x(t) = a + bt + / h(t - r)flf(r)n(r) dr. 

J— 00 

We define b, K(t, z) and h(l — t) as in Case 1 and again require that' 
in the absence of noise, 6 equal b, or that the constraints (15) hold. 
These assumptions lead to the following integral equation: 

/ dz'Kit, z')${t - z,t - z) = X + pz, 

Jo 

where 

*(t - z,t- z') = [ h(ij)h(y + z' - z)g\t - z - y) dy. 
Jo 

Instead of dealing directly with the integral equation, we first rewrite 
x(t) as 



x(t) = f h{l - t)[0 + yr + g(r)n(r)] dr, (19) 

J— 00 



where 







7 = 


a b 
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Jo 
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Jo 
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from the assumed form of T[h(t)], and that 
/ h(t - t)(/3 + yr) dr -= 

J— 00 

»oo /•" r 00 

d / h(x) dx + 7* / /i(.r) f/.c - t / xh(x) dx = a + hi 

Jo Jo Jo 

when the values of /3 and 7 are substituted from (19). 

A function y(t) is now defined as 

y(t) - /3 + 7' + <?(*)"(*) 

and a function A'(/, 2) is defined as the convolution of K(t, z) and h(z), 
so that (14) is replaced by 



b = I K(l, t - t)ij(t) dr. 

J— CO 



(20) 



The function y(t) is similar to x(t) as defined in (2). Therefore, K{t, z) 
may be found in a manner analogous to that used in the case in which 
the noise was not pre-filtered. Recalling that K(l, z) is required to 
vanish outside (0, T), we may define K{t, z) more exactly by 

K(t, z) = f dvK{t, v)h{z - v); < z < T 

Jo 

= f dvK(t, v)h(z - v); z> T 
Jo 

= 0; z < 0, 

or the equivalent: 

" tV '" d' 

$„!,*«.«>-£«,£«<,.); o<,<t (21) 

= 0; z > T, z < 0. 

Since R(t, z) may be discontinuous at 2 = and T, these points have 
been excluded. While /v(/, z) is constrained by (15), it follows that 
K(t, z) must satisfy 

I K(t, z) dz = 

Jo 

(22) 
f zK{t,z)dz = -^. 

Jo «0 
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The variance to be minimized is now 

E(b - b)' 2 = j K\t, z)g\t - z) dz. 
From (21), 

R(t, z) = £ a /€ fjS ; z > T, (23) 



where the Sj are presently unknown but will be selected to minimize 
E{b — b) and the f ; are known in terms of the <xj . A variation of 

[ [K\t, z)g\t - z) + XK(t, z) + nzK{t, z)] dz 

Jo 

gives 

2K((, z)g\t - z) + X + fiz = 0; < z < T (24) 

and a set of n equations defining the Sj : 

n — 1 ,.00 -oo 

2 Z sj I e lti ™'g\t - z) dz + / (X + ^/jf tfe = 0, (25) 

where 

k = 0, 1, .-.,n - 1. 

Equations (23) and (24), together with the constraints (22), define 
K(t, z) completely. The function K(t, z) may now be found from (21). 
From the discontinuities in R(t, z) at z = and T, delta functions 
of order as high as n — m — 1 may be expected in K(t, z) at z = and 
7\ For example, if A(.r) equals e~ ax , (21) becomes 
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which defines s and therefore K(t, z) for z > T. Except for evaluating 
X and m from (15), the solution is 



K(l, z) = U 



K-L \ ;,■■.) -8(z-T)+a+~ 



X + \xz 



X . 

aH(t, T) 



+ 






where 



//(/, T) = e 2ar /" e -*tfG - 2) tfe. 
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